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ABSTRACT 

The study of multiple extrasolar planetary systems has the opportunity to obtain 
constraints for the planetary masses and orbital inclinations via the detection of mu- 
tual perturbations. The analysis of precise radial velocity measurements might reveal 
these planet-planet interactions and yields a more accurate view of such planetary sys- 
tems. Like in the generic data modelling problems, a fit to radial velocity data series 
has a set of unknown parameters of which parametric derivatives have to be known by 
both the regression methods and the estimations for the uncertainties. In this paper 
an algorithm is described that aids the computation of such derivatives in case of 
when planetary perturbations are not neglected. The application of the algorithm is 
demonstrated on the planetary systems of HD 73526, HD 128311 and HD 155358. In 
addition to the functions related to radial velocity analysis, the actual implementation 
of the algorithm contains functions that computes spatial coordinates, velocities and 
barycentric coordinates for each planet. These functions aid the joint analysis of multi- 
ple transiting planetary systems, transit timing and/or duration variations or systems 
where the proper motion of the host star is also measured involving high precision 
astrometry. The practical implementation related to the above mentioned problems 
features functions that make these kind of investigations rather simple and effective. 

Key words: Celestial mechanics - Methods: Analytical, Numerical - Methods: iV- 
body simulations - Techniques: radial velocities 



1 INTRODUCTION 

As of this writing, 36 multiple planetary systems are known 
around main sequence stars. Most of these systems bear two 
detected planets, 8 of them have 3 and 2 of them have 4 
planets while the star 5 5 Cnc has 5 compan ions . With the 
exception of HR 8799 (|Marois et al.ll2008h . all of the de- 
tections are based on or confirmed by the measurements of 
the radial velocity (RV) variations of the host stars 2 . The 
planet HAT-P-13b l|Bakos et al.ll2009r ) also transits its host 
star. The detection of this planet was based on transit pho- 
tometry while the second companion in this system has been 
revealed by radial velocity measurements. In general, analy- 
sis of RV variations constrains the mass (m) of planets by a 
lower limit. Namely, only the quantity msini is determined 
by RV data where i is the orbital inclination (relative to 
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See e.g. http://exoplanet.eu for an up-to-date list. 
2 The planetary system around HR 8799 with 3 confirmed planets 
has been detected by direct imaging. 



the tangential plane of sky). With the exception of transit- 
ing planets and systems were the spatial motion of the host 
star is detected via astrometry, there is no direct evidence 
for the actual value of the orbital inclination (and therefore 
the mass of the planet). In case of transiting systems, incli- 
nations are constrained by m easuring the im pact parameter 
from the light curves (see e.g. IPal et al.ll201CT ). while astrom- 
etry yields not only t he inclination but the o r ientation of the 
orbital plane as well (Bean fc Seifahrtll20091 : iBenedict et al.l 



I201CJ : iMcArthur et alj|2010h . The planetary system around 
GJ 876 is the only known one whe re mutual inclination is 
also detected with a 2-<r confidence (|Bean fc Seifa hrt 2009). 
A great advantage of multiple planetary systems is the pos- 
sibility of detecting mutual perturbations via the deflection 
of RV values from the pure ly Keplerian solution (see e.g. 
lLaughlin fc Chambers! 1200 ll ). Therefore, precise analysis of 
accurate RV series may yield to an acceptable constraint for 
both the inclination and the planetary masses. Addition- 
ally, planet-planet interactions depend on the mutual incli- 
nations, thus more complex models (such as non-coplanar 
orbits) for the whole planetary systems can be investigated. 
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Like in the majority of data modelling problems, RV 
variations (of the host star) in single or multiple plane- 
tary systems are modelled with a function of a few exter- 
nal (unknown) parameters. These parameters include the 
orbital elements and the masses of the planets as well as 
the barycentric velocity of the host star. For most of the 
regression methods involved in data modelling (se e e.g. the 
Levenberg-Marquard algorithm. IPress et al.lll992f ). and for 
the analytic estimation of the covariances, uncer tainties and 
correlation s of the model parameters (see e.g. iFinnl Il992l ; 
|Palll2009al ). the partial derivatives of the model functions 
(with respect to the model parameters) have to be known in 
advance. The simplest way of RV curve modelling does not 
take into account the mutual interactions between the plan- 
ets and characterizes the observe d RV variations as a sum 
of independent Keplerian models. IWright fc Howard! (2009) 
describes an algorithm detailing the efficient computation 
of the parametric derivatives of RV model functions where 
the mutual planetary perturbations are neglected. The main 
objective of this paper is to present an algorithm that cal- 
culates parametric derivatives when the planet-planet inter- 
actions are also taken into account. 

The structure of the paper is as follows. In the next sec- 
tion, we describe the mathematical tools used to construct 
the algorithm itself (including the discussion of optimal or- 
bital parameterization, and the numerical integration). In 
Section [3] the practical implementation is detailed while in 
Section [4] we demonstrate the usage of the algorithm for 
three specific multiple planetary systems. The results are 
summarized in the last section. 



2 ORBITAL PARAMETERIZATION AND 
NUMERICAL INTEGRATION 

In this section we summarize the conventions involved in the 
orbital parameterizations (as used throughout this paper), 
the algorithms used for numerical integrations and the meth- 
ods applied in the calculations of parametric derivatives of 
the radial velocity model functions. 



2.1 Spectroscopic and Keplerian orbital elements 

The measured radial velocity variations of the host star de- 
termine the spectroscopic orbital elements of the planetary 
companion(s). In the case of a system with a single planet, 
there are 6 of such parameters: the period, P, the time of 
periastron passage, T p , eccentricity, e, argument of perias- 
tron, w, the semi-amplitude of the RV variations, K, and 
the zero-point (or the mean) radial velocity of the central 
body, 7. The above set of orbital parameters is the most 
widely used in the literature, however, it is not the best 
choice for our purposes because of the following reasons. 
First, for nearly circular orbits the argument of pericenter 
is not so well constrained and for exactly circular orbits it 
cannot be defined at all. Second, the time of pericenter pas- 
sage is also purely constrained for orbits with small eccen- 
tricities and not denned for e = 0. In order to avoid such 
ambiguities, one can use the Lagrangian orbital elements 
k — e cos ui and h = e sin ui for the parameterization of the 
shape of the orbit and use one of the following quantities 
instead of T p : the mean longitude A = \(Eo) or the orbital 



longitude tp = ip(Eo) for a certain fixed epoch Eo, or the 
moments T\ or T Vo when the mean longitude or the orbital 
longitude have a certain fixed value of Ao or ipo, respectively. 
All of the above four quantities are well-defined for circular 
and nearly circular orbits. The element T VQ is widely used 
in the case of transiting planets since due to the definition 
of the alignment of the reference frame, transits occur at 
<Po = 7r/2. In the case of multiple planetary systems where 
the mutual interactions are not negligible, the usage T\ Q 
or T vo is not the best choice since in practice these imply 
different epochs for the distinct planets. Throughout this 
paper, we use the mean longitude A as the primary orbital 
element to characterize the phase of the orbital motion. The 
conversion between T p and A is rather simple, namely 

2tt 



A 



(E - T p ) + u. 



(1) 



Since the mean longitude at an arbitrary moment t is 
A(i) = n(t — So) + A, i.e. it is a linear function of the mean 
motion n — (2ir)/P and the A, we prefer n to P. There is an 
additional benefit using n for characterizing orbital periods: 
simple error propagation estimations yield that for a given 
RV semi-amplitude, the uncertainty of n is independent for n 
itself. In other words, in a multiple planetary systems where 
planets have masses with the same magnitude, one can ex- 
pect that the obtained uncertainties are also roughly the 
same. This is also confirmed by the demonstration analysis 
presented here (see Sec. 14.1(1 . 

In order to interpret the spectroscopic orbital elements, 
these should be converted into Keplerian parameters. For 
any type of orbits, the orbital eccentricity, argument of pe- 
riastron, and the mean longitude are interpreted in the same 
way. In the case of planar orbits, semimajor axis a of the or- 
bit and the mass of the planet m are derived from the mean 
motion n, and the normalized RV semi-amplitude K, = KJ, 
where J = y/1 — e 2 = \/l — k 2 — h? . Let us denote the mass 
of the central star by M . Using Kepler's Third Law and the 
barycentric velocity of the central body, namely 

aV = G(M + rn) 



and 

K. — an- 



M+m 

we obtain for Gm and a: 



Gm 



(GM)F 



£ 3 



\GMri, 
2/3 {GM + Gm} 1 ' 3 



(2) 
(3) 

(4) 
(5) 



-2/3 



GM 



( K 3 



\GMn 



1/3 



Here F m (a) denotes the solution of the equation 
x 3 

(1 + %y 

If a ^ 0, the above equation always has a unique solution. 
Moreover, if a > 0, the F m (a) function behaves analytically 
and its derivative is 
dF m (a) [l + F m (a)]F m (a) 



da 



a[3 + F m (a)] ' (7) 

The well-known proportionality K, oc Gm is a direct conse- 
quence of F m (a) ~ a 1 ' 3 for a«l. We note that in the case 
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of spatial orbits when the orbital inclination i differs from 
90°, the gravitational parameter and the semimajor axis are 
calculated as above, but the normalized RV semi- amplitude 
is /C = i^J(sini) _1 . 

Once the Keplerian orbital elements are derived, the 
planar coordinates and velocity vector components of the 
planet can be calculated as 



+ 



+h 

l + J\-k 



—s 

+c 



+h 

l+J\-k 



(8) 



(9) 



Here q — ecosE, p — esinE 1 , c = cos(A + p) and s — 
sin(A + p), while E den otes the eccentric anomaly. As it 
is shown by IPall ([2009a), the quantities q, p, c and s are 
analytic functions of the mean longitude and the Lagrangian 
orbital elements (k,h). The observed radial velocity of the 
central body is then 



K (1) = - 



M + m 



(oxX + Oyy), (10) 



shows the direction of the 



where the unity vector (o x 
observer. Due to the historical definition of the reference 
frames used in the astrophysics of binary stars and extraso- 
lar planets, this vector is fixed to be (o x , V ) — (0, 1). 



2.2 Lie-integration 

In this short section we summarize the properties of the nu- 
merical integration of ordinary differential equations named 
Lie-integration. The main featur e of the numerical inte- 
gration based on the Lie-s e ries (IGrobner fc Knapp 119671 : 
iHanslmeier fc Dvorakl 1 19841 : lEggl fc Dvorakl I2010T ) is that 
the solution of the differential equation 



/i(x). 



(11) 



is approximated by its Taylor series up to a finite order. The 
coefficients for this power series expansion are generated by 
the Lie-operator 



U := J2 f^ 



(12) 



(where Di = -^-) with which the solution of equation (|lip 
can be written as 

x(i + At) = exp (At • L Q ) x(t) 

where 



(13) 



exp( 



Af rfc 
fc! 



At ft 



(a^o) = £^ = £^t 5>a 



(14) 



The method of Lie- integration is the finite approximation of 
the sum in the right-hand side of equation (|14p . up to the 
order of M. Thus, the solution after At time is approximated 
by 



At fe 



(* + Ai H E ^r L o *w = E it ( L « x(t) ) • (15) 



fc=0 / fc=0 

fc. 



Supposing the coefficients L x.(t) are computed, the numer- 
ical integration itself is straightforward. In practice, the val- 
ues of these coefficients are computed involving recurrence 



relations, i.e. for n ^ 0, the Lq +1 x(£) terms are evaluated 
using the previously calculated L§x(£) (0 ^ k ^ n) coeffi- 
cients. For each particular problem, the recurrence relations 
must be derived properly. For the general JV-body problem , 
these relations are presented in lHanslmeier fc Dvorakl (|l984r ) 
while the relations in the cases when t he reference frame is 
fixed to one of the bodies are shown in lPal fc Sulil (J2007J). 

The computational cost (CPU time) required by the 
calculation of the LQx(t) coefficients is de finitely larger tha n 
the time of evaluating equation (|15p (see lPal fc Siilil 120071 ). 
Thus, a great advantage of the Lie-integration is that the 
stepsize At can be altered after the coefficients are eval- 
uated without much of additional cost, and therefore an 
effective integration method can be implemented with an 
adaptive stepsize. Moreover, there is an availability of an 
alternate way of adaptive integration, namely if the stepsize 
At is kept fixed, the integration order M can also dynami- 
cally _be_JiiCTeased^ntilwe reach the desired precision (see 
also lEggl fc Dvorakll2010r h More details about the practical 
implementation are found in Sec l3.ll 



2.3 Motion in the reference frame of one of the 
bodies 

In several applications, such as in the description of a plan- 
etary system or in perturbation theory, the equations of 
motion are transformed into a reference frame whose ori- 
gin coincides with one of the bodies. Practically, this is the 
body with the largest mass, i.e. in a planetary system it is 
the host star. Let us define the central body as the body 
with the index of i = 0. Altogether we have 1 + N bodies, 
where the other ones are indexed by i = 1, . . . , N. Let us 
use the relative (non-inertial) coordinates Tim and velocities 
Wi m (where the second index m refers to the spatial dimen- 
sion, i.e. m = 1 is for the x coordinate, m = 2 is the y 
coordinate and in non-planar problems, m = 3 refers to the 
z coordinate). The equations of motion in a compact form 
are 



Tim — ^jm , 

w im — —G(M + mi)4>ir ir , 

N 



G ^ rUj [<f>ijAijm + </}] r j 
3 = 1,3^1 



(16) 



(17) 



where Aij m is the mth component of the vector pointing 
from the body j to body i, 



^ijm - — Tim Tjm> 



(18) 



Pi and pij denotes the distances from the central body and 
the mutual distances, respectively: 



pij : — / / j AjjmAjjm, 



poi = pm = i^2n m n m , (19) 

(20) 
and (pi and cf>ij are defined as the reciprocal cubic distances, 

<ki ■■= pTj 3 , (2i) 

& := P~ 3 - (22) 
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Note that the quantities pi and pij, like so <pi and 4>ij are 
distinguished only by the number of their indices. With- 
out going into details, we present the recurrence relations 
of the Lie-derivatives, including the linearized variables in 
Appendix [Al 

Using the notations introduced above, the observed ra- 
dial velocity is computed as 



V T 



(N) 



E 



/ N v 

J2 miWim 



N 

£ 



N 



(23) 



where the unity vector o = o m — (°x,o y ) or o m = 
(o x ,o y ,o z ) defines the direction of the observer (see also 
equation [TO)) . 



2.4 Parametric derivatives 

In the analysis of radial velocity variations in a multiple 
planetary system, the parametric derivatives of the RV curve 
have to be computed at a certain moment t with respect to 
the initial orbital elements (at t = Eq). Let us define an 
arbitrary quantity Q which depends only on the solution of 
the ordinary differential equation (|11[) . Q = Q(x(i)) = Q(x). 
The parametric derivatives of the quantity Q with respect 



to the initial conditions x 



o 



can be computed in the 



following way. First, let us write the linearized equations of 
equation (|11[) as 



6=6 



dx m 



(24) 



(Here and in the following we use the implicit summation 
notation wherever it is unambiguous.) The variables £& de- 
note the so-called linearized variables for one particular ini- 
tial condition. Second, due to the linear property, with the 
solution of the full linearized equations 



Zfh — Zf. 



g/Ux) 

dx m 



(25) 



one can compute the solution of equation (|24p for any arbi- 
trary initial conditions £°, namely: 



&(t) = -Z M (i)6° 

if the respective initial conditions of equation (|25|) are 



Zi , 



= 6th = 



if l = k, 
if I ^ k. 



(26) 



(27) 



Finally, without going into the details, it can be shown that 
the partial derivatives of Q(t) with respect to the initial 
conditions x° = xL „ is 



dQ(t) 

dx° f 



Zih{t)— — . 

OXk 



(28) 



where Zik (t) represents the solution of equation (|25p at the 
instance t. If the initial conditions x° are defined with an 
alternative parameterization, i.e. x° = x°(x), the paramet- 
ric derivatives of the quantity Q with respect to the x are 
calculated involving the chain rule, namely 



dQ 

dx e 



dx° e dQ 



(29) 



In the analysis of RV data series, x represents the set of 
spectroscopic orbital elements - including the (normalized) 
RV semi-amplitude, x = (JCi,m, Xi, ki,hi) -, x° represents 
the spatial coordinates, velocities 3 and the gravitational pa- 
rameters of the planets, while Q = V T , the observed radial 

dx° 

velocity. In practice, the partial derivatives -^-A- can be com- 
puted using the formulae presented in Appendix|B] while the 
computation of the terms -*-£- is relatively simple since in 
the equation for the radial velocity (see equation [TD] or later 
in Sec. 12.31 equation [23} is a rational expression of two func- 
tions in that are linear with respect to both the coordinates 
and masses. 



2.5 Linearized equations 



As we have seen before (Sec. 12. 4f) , linearized equations have 
to be solved in order to calculate the partial derivatives of 
an arbitrary quantity (that depends on the solution of the 
original differential equation) w ith respect t o the initial con- 
ditions. As it has been shown in lPal fc Siilil i|2007l ). using the 
same notations as above the Lie-derivatives of the partial 
linearized variables £j (see also the previous subsection) can 
be written as 



L £h — £,mD m L Xk — (,mD m L Xk- 



(30) 



Obviously, this formula can be applied to obtain the solution 
for the full linearized form (see equation I25[l : 

L Ztk = Zt m D m LQXh. (31) 

Thus, the solution for equation (|31[) has to be substituted 
into equation (|28[) or (|29p in order to obtain the partial 
derivatives of arbitrary quantities with respect to the initial 
conditions. The complete set of recurrence relations for the 
linearized problem is found in Appendix |A"1 



3 IMPLEMENTATION 

The algorithm presented in Section [2] has been imple- 
mented as an add-on module for the regression analysis 

anol d ata m odelling program If it 4 (described briefly 

in HU l2009bl ) in a form of an ANSI C code. Since the 
application program interface (API) of If it for this 
kind of dynamically loaded libraries is rather simple, the 
source code module can easily be modified for arbitrary 
purposes, such as inclusion for other kind of C programs 
or another languages or programming environments 
that support linking of C modules (e.g., FORTRAN or 
IDL). The source code is available from the web ad- 
dress http: //szof i . elte.hu/~apal/utils/astro/nbrv. 
A standalone implementation of the Lie- 
integrator code is also available from the address 
http: //szof i . elte.hu/~apal/utils/astro/lieint, 
with the same algorithmical features and with an easy 
user interface for numerical integration and simple stability 
investigations. 



3 See also Sec. l2.3l for further details about the notations used in 
the description of multiple planetary systems. 

4 This program is available as a part of the libpsn package, see 
http : //szof i . elte . hu/~apal/utils/libpsn/. 
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Table 1 . List of additional modes as implemented in the generic functions nbrv_2g_N ( ) and nbrv_3g JJ ( ) . These functions have 4 or 5 
additional parameters (in the respective cases of 2 and 3 dimensional variants) comparing to the pure RV functions (nbrv_2d_N ( ) and 
nbrv_3d_N ) . The first additional parameter is the "mode flag", F, an integer between and 3. The second parameter is the body 
index, k, a non-negative integer less than or equal to N. The other 2 or 3 parameters are the components of the o m vector. All of the 
values computed by these functions are projected coordinates or velocities: the spatial vectors are multiplied by the o m components, 
yielding a scalar product. Thus, in this table such derived coordinates and velocities are referred as "projected coordinates" or "projected 
velocities" . 



Mode Body 
(F) index (fc) 



Result 



Interpretation and typical usage 



V W 

' r . m u n 



VimO m Projected velocity of the barycenter with respect to the central body. In case of planetary 

systems, interpreted as the radial velocity of the host star where the line-of-sight is defined 
by the o m vector. If o = (0, 1) or o = (0, 0, 1), the results are equivalent with the results of 
nbrv_2d_N() and nbrv_3d_N ( ) . 

(N k) 

1 ^ k ^ N V r ,m' o m Independent components of the projected barycentric velocity. Although only the joint effect 
of all of the planets in the planetary system can be measured by radial velocity variations, 
the contribution of each planet to the final RV curve can be analyzed by this way. Due to 
the mutual perturbations, these velocity components are not strictly periodic and cannot be 
described only by the orbital elements of the respective planet. 

8; im 'o m Projected coordinates of the barycenter with respect to the central body. In case of planetary 
systems, these coordinates describes the wobbling of the host star, as it might be detected 
by precise astrometric measurements. 

(N k) 

1 ^ k ^ N -Br,m o m Independent components of the projected barycentric coordinates. If the complementary 
inclination for a particular planet is close to zero, the planet transits the host star, yielding a 
small flux decrease that can be measured. In case of suchtransiting planets, these coordinates 
determines the magnitude of the transit timing variations due to light-time effects. 

1 ^ k ^ N w^ m o m Projected spatial velocity of the body k. For planets with a nearly edge-on orbit, the tangen- 

tial acceleration of transiting planets is negligible at the time of the transits. These velocities 
well constrain the duration of these transits. 

1 ^ k ^ N r^mO-m Projected spatial coordinates of the body k. These coordinates constrain the shape of these 

possible transit light curves as well as the precise timings of the transits if an orbit is nearly 
edge-on. 



Table 2. Spectroscopic orbital elements for the planetary systems HD 73526, HD 128311, and HD 155358. These orbital elements have 
been derived from the data available in the literature (see text for further references) and have been used as an initial condition in the 
fits discussed in this paper. The last column shows the number of available radial velocity data points (the same as involved in the fits). 



System 


E (BJD) 


Ah/M Q 


Ki 


sin i 


(m/s) ': 


%i = 2n/P l (1/d) 


A, 


(rad) 


k L 


= ei cos u>i 


hi 


= ei sintJi 


AW 


HD73526 


2,452,500 


1.08 






70.0 
61.4 


0.03360 
0.01620 




3.902 
4.150 




-0.402 
-0.480 




+0.040 
-0.080 


31 


HD128311 


2,452,500 


0.84 






64.6 
75.1 


0.01370 
0.00677 




1.896 
1.500 




-0.090 
-0.160 




+0.233 
-0.058 


75 


HD155358 


2,453,500 


0.87 






34.6 
14.1 


0.03222 
0.01185 




0.894 
0.249 




-0.106 
+0.027 




+0.035 
-0.174 


71 



Table 3. Best-fit spectroscopic orbital elements and their 1-ct uncertainties for the planetary systems HD 73526, HD 128311, and 
HD 155358 derived by the MCMC algorithm under the assumption of coplanar orbits. The median for each probability distribution is 
treated as a best-fit value. See text for further details. 



System 



/Qsini (m/s) rii = 2ir/Pi (1/d) 



Xi (rad) ki 



1 (m/s) 



HD73526 



65.9 + 3.8 
61.8 + 0.5 



0.03359 + 0.00017 
0.01629" 1 



-,+0.00018 
-0.00015 



3.855 + 0.046 
4.108 + 0.076 



-0.404 ± 0.045 
-0.503 



0.044 
0.041 



+0.069+g;8£ 

-0.046 ± 0.046 



-33.6 + 2.4 0.82i^ 



HD128311 



49.1 



+9.2 
-5.9 



0.01364 + 0.00012 



1 5OQ+0219 



73.8 + 3.1 0.00663 + 0.00010 1.396 + 0.046 



+0.0351°;^ 

in 142+0-093 
+ U14/ -0.169 



+0 337+°- 0n 

-T-U..M ( _g 07g 

+0-063ig;g83 



0.9 + 2.1 0.8toJ 



HD155358 



31.18 + 0.34 
13.66+1?° 



0.03225 



+0.00024 
0.00016 



01203+ 000029 
u - ulzuo -0. 00033 



0.867 ± 0.064 
161+ 0147 

U.1D1_ 169 



-0.136 + 0.039 
-0-065tnfs 



+0 041+ 0043 
-l-u.u4i_ 047 

-0 138+ 0116 

U.13O_ 179 



10.1 + 1.0 
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In practice, this add-on module (named nbrv . so on 
most of the UNIX systems or nbrv.dylib on OS/X) reg- 
isters the functions named nbrv_2d_N() and nbrv_3d_N() 
where N = TV is the number of planets in the planetary 
systems 5 . These functions have 1 + 5TV + 1 or 1 + 7TV + 1 pa- 
rameters, for the respective cases for the planar (nbrv_2d_M) 
and spatial (nbrv_3d_N) problems. The first parameter is the 
central mass M. (in Solar units), the following TV x 5 or 
TV x 7 parameters are the spectroscopic orbital parameters 
K.i (in the units of m/s), n* (in the units of d _1 ), \i (in 
radians), k% and hi. Furthermore, the spatial functions have 
two additional parameters: the complementary orbital incli- 
nation, i = 90° — i and the argument of ascending node, Q, 
both angles are measured in radians 6 . All of these orbital 
elements are defined for a certain epoch of Eq (BJD). The 
last parameter is the time At elapsed from the epoch Eq, 
i.e. At — t — Eq. Here i is the index for the actual planet, 
1 < i sC TV. 



3.1 Adaptive integration 

As mentioned earlier, one of the advantages of the Lie- 
integration is the possibility of the implementation of a two- 
way adaptation, by varying both the stepsize and the or- 
der of the integration. Since the summation of the Taylor- 
coefficients in equation (|15|) does not need so much comput- 
ing time (compared to the evaluation of these coefficients), 
the stepsize of the integration can easily be altered in or- 
der to reach the desired precision. In practice, the adaptive 
integration is implemented as follows. First, let us define a 
minimal and maximal order of M m i n and M max . It is easy to 
see that an initial, nearly optimal stepsize for the integra- 
tion is Aio oc n„a X , where n max is the maximum of the mean 
motions appearing in the planetary system. In the case of 
circular orbits, Ato = 0.8 n max is a good choice for M w 20 
and for a relative precision of 5 = 2 ■ 10 -16 . If we allow 
a minimal and maximal order for the Lie-integration (M m i n 
and M max , respectively), the adaptive control of the stepsize 
and integration order is done as: 

1. The terms LqX in equation (|15[1 are evaluated using 
the appropriate recurrence relations and the summation is 
performed with a fixed value of At. 

2. If the desired precision 8 is reached before the order of 
M — Mmi n , At is multiplied by the factor M max /M m i n and 
re-compute the sum in equation (I15|l (with the additional 
evaluation of the necessary terms where k > M). This step 
is repeated until M < M m i n . 

3. If the desired precision 5 cannot be reached until the 
order of M — M max , divide At by the factor of M max /Af m i n 
and re-compute the sum in equation (|15[) . This step is re- 
peated until M max < M. 

4. If the desired precision is reached between the orders of 
M m i n and M max , accept the value of At and proceed with 
the next step of the integration. 



We have to note that this kind of two-way adaptive inte- 
gration assures that we definitely obtain the desired preci- 
sion level without the loss of computing time. In the case 
of more common integrators (like Runge-Kutta or Bulirsch- 
Stoer methods) , the estimation of the accuracy is based on 
heuristics and it is not checked by these algorithms that the 
expected precision is really obtained. If it turns out that the 
stepsize is too large (or other parameter of the integration 
should be changed), then a total re-computation is needed 
for these algorithms and we definitely lose the computing 
time spent on the previous evaluations. 

In the actual implementation of the nbrv module, 8 has 
been chosen by default to be the precision level of the IEEE 
double precision (64 bit) floating point number representa- 
tion, that is 8 — 2-10~ 10 . Although it is an extreme precision 
compared to the implied and required precision level for the 
problem, this precision implies that the whole set of func- 
tions implemented in the module can be treated as analytic 
functions without any side-effects. Additionally, this preci- 
sion level ensures that there would not be any systematic 
distortions by varying the samples on the domain of inves- 
tigations. 

3.2 Properties 

In the following, we present some properties for these 
functions. For simplicity, let us denote by S; the set of 
the spectroscopic orbital elements (ni,Xi,ki,hi) and define 
Vi (•) = nbrv_N(.). Since interaction between the planets 
is relatively small (comparing to the gravitational force of 
the central star), 

V} N) (GM , JCi , Si , . . . , IC N , S N , At) « 



^y r (1) (GM,/C„S„At). (32) 



It is easy to show that for fixed values of Id, the effect of 
mutual interactions decreases as the central mass increases, 
namely 

lim V I (N) (GM,fCi,S 1 ,...,K, N ,S N ,At) = 

GM.— too 



^2v r (1) {GM,}Ci,Si,At). (33) 



Similarly, scaling of the normalized RV semi-amplitudes 
yields the same kind of equation: 



lim (sin i) V}" > GM , -^- , Si 



sin %~ >oo 



K.i 



/C/> 



= ^K (1) (G7W,/C l ,S l ,At). (34) 



5 In the current implementation N ^ 8, however, the source code 
can easily be modified to increase the maximum number of plan- 
ets. 

6 The complementary angle of the inclination is used for simplic- 
ity: the planar functions yields the same values as the spatial ones 
if these two additional parameters are set to zero. 



Although | sin i\ ^ 1 for real orbits, either the above equation 
or equation (|33|) can be used to for mally decrease the level 
of in teraction between the planets IjLaughlin fc Chambers! 
l2001f l. Additionally, the fit of independent Keplerian orbits 
involving the right-hand side of equation (|32l) yields good 
initial conditions for the real problem. 
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Additionally, using equations © and (|l(jp , V r can also 
be written as 



V T w (GM,)Ci,Si,At) 



1-q 



cos(A + p) 



kiq 

l + J 



(35) 



where Si = (m, Ai, ki,hi), A — n\At + Ai, p — p(A, ki,hi), 
q — q(A, ki,hi) and J = \J\ — kf — h\. Obviously, V T does 
not depend on GA4 , therefore this argument is only a formal 
one. 

As it is known from the literature of binary stars and 
hierarchical stellar systems, if JV = 1, variations in the ar- 
gument of the ascending node have no observable effect and 
like so, for two planets, radial velocity variations depend 
only on the difference D = Q2 — Sli (that also determines 
the mutual inclination of the two orbit). In general, for all 
N ">■ 1, we can state 



E 



dV T 



(JV) 



dQ k 



(36) 



that is also equivalent for the previously mentioned special 
cases for N = 1 and N = 2. One should keep in mind 
these properties while utilizing the nbrv_3d_N() functions 
for purely radial velocity data. 



3.3 Generic functions 

Throughout this paper we were focusing on the analysis of 
radial velocity variations of host stars in (multiple) plane- 
tary systems. However, the results of the numerical jV-body 
integrations can be exploited in another aspects of binary 
and hierarchical stellar systems and extrasolar planet studies 
- including eclipsing binaries, transiting planets and plane- 
tary systems were astrometric information is also available. 
For instance, using the spatial coordinates rum, one can esti- 
mate the moments of eclipses and transits and characterize 
the shape of the light curves. Similarly, the velocities Wkm 
can directly be used to calculate the durations of the eclipses 
and transits, since at the time of these events, the tangen- 
tial acceleration of the transiting body is nearly zero, so the 
tangential velocity computed from Wkm is a rather good ap- 
proximation for the reciprocal duration 7 Obviously, using 
the capabilities of the program If it, the knowledge of the 
velocities and coordinates can be exploited in complex stud- 
ies were simultaneous fits are performed on, for inst ance, RV 
and astrometric data series IjMcArthur et alj|20ich . 

First, let us write the independent components of the 
barycentric coordinate and velocity in the form of 



B 



(N,k) 



v T 



(JV.fc) 



N 

i=l 

rUkWkm 

N 

i=l 



and 



(37) 
(38) 



It can be shown that the projection of the barycentric coor- 
dinates to the line-of-sight is proportional to the light-time 

7 In practice, total transit duration must be computed by taking 
into account the impact parameter that is derived from the orbital 
inclination and the alignment of the orbital ellipse. Therefore, the 
observed duration depends on the rj. m coordinates as well. 



effects, thus in the analysis of timing variations in eclipsing 
and/or transiting systems, these corrections should also be 
taken into account. With the terms defined in equation (1381) . 
the observed radial velocity of the host star can be written 



v; 



(JV) 



= £ 5> 



r(N,k) 

" r.m 



(39) 



The independent components V r \ m ' show the influence of 
each body on the final observed radial velocity variations. 
However, due to the mutual interactions, these functions are 
not strictly periodic. 

In order to analyze the effects discussed above in mul- 
tiple stellar and planetary systems the module nbrv . so im- 
plements the functions nbrv_2g_N() and nbrv_3g_N ( ) (where 
N = N is the number of planets). These functions have 
4 or 5 additional parameters 8 comparing to the functions 
nbrv_2d_N() and nbrv_3dJJ(). The first additional param- 
eter is a mode flag, that is used as a selector between the 
V T , B T , w and r quantities. The second additional parameter 
is the body index k while the last 2 or 3 numbers represent 
the (o x , o y ) or (o x , o y ,o z ) vector. The value computed by the 
nbrv_2g_N() and nbrv_3g_N() functions is a scalar product 



of the quantities r kr . 



U'kn 



B, 



(JV,fc) 



and V, 



(JV,fc) 



and the o„ 



vector. The comprehensive list of the modes implemented in 
the functions nbrv_2g_N() and nbrv_3g_N() can be found in 
Tabled] 



4 APPLICATIONS 

In this section we describe two possible applications of the 
algorithms presented in this paper. First, we show some ex- 
amples related to the problem of orbital characterization. 
In the second part, we describe in brief how optimal obser- 
vation scheduling can be performed by employing the If it 
code and this implementation discussed above. 



4.1 Orbital fits 

As a demonstration and application of the algorithm 
described in Sec. [2l we analyze the radial veloc- 
ity data fo r the multip le (double) planeta ry systems 



lty data to r tne multip le (doublej planeta ry s ystems 

HP 73526 dTinnev et all I2006T ). HD 1283 11 dVogt et al 



120051 ; IWright et al.l l2009h and HD 155358 l|Cochran et al 



[2007J) in the form as it is available from the literature. The 
objective of our test presented here is to constrain the or- 
bital inclination assuming a planar model for these plane- 
tary systems. Coplanar orbits are expected from modelling 
l|Goldreich. Lithwick fc Sarill2004h and also confirmed by ob- 
serva tions (Solar System or GJ 876, see iBean fc Seifahrd 
120091 ). We note here that complex dynamical and stabil- 
ity investigations for multiple planetary systems can also 
be used to constrain or refine orbital ele ments (see e.g. 
iFerraz-Mel lo. Michtchcnko, & Bcauge 2005), imply alterna- 
tive planetary configurations (]Gozdziewski fc Konackill 20061 ) 
or rule out small in clinations (in parallel wi th direct RV data 
modelling, see e.g. lLaskar fc Correiall2009l ). 

We employed the method of Markov Chain Monte-Carlo 



In the cases of 2 and 3 dimensional variants, respectively. 
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Figure 1. Probability distributions of the sini orbital element for the planetary systems HD 73526 (left panel), HD 128311 (middle 
panel) and HD 155358 (right panel) derived from the radial velocity data available from the literature. 



(MCMC, see e.g. lFordll2005t) in order to derive the best-fit 
orbital elements for these interacting planetary systems. The 
initial conditions (orbital elements) were based on the data 
available in the literature and summarized in Table [2] This 
table shows the orbital elements for a certain epoch (-Eo) 
while the eccentricity, argument of pericenter and the time of 
pericenter passage have been converted to the spectroscopic 
orbital elements discussed in Sec. 12.11 In all of the fits, the 
model function 



RV 



7+(sirn)V; <2) ( 



(/M,4^,S!,4^,S 2 ,t- 
sm i sm i 



E ) (40) 



has been utilized. In the MCMC runs, the values of sini 
have been forced to be between 0.2 and 1 (note the prob- 
ability that sini of a randomly oriented orbit is less than 
0.2 is roughly 2%). The results of the fits are displayed in 
Tabled] The a posteriori distributions of the sini values are 
shown in FigQ] As it is clear from the plots, in the case of 
HD 73526, the orbital inclination is well defined purely by 
the RV data, although the upper limit for sini is not con- 
strained. Namely, 0.68 < sini for this system within 1-a, 
that is equivalent with 42° < i. For the other two planetary 
systems, the available RV data do not provide a significant 
constraint for the inclination. 

An advantage of the knowledge of the partial derivatives 
of the model functions that the uncertainties of the fit pa- 
rameters can be estimated analytically involving t he method 
of Fisher matrix analysis (JFinnlll992i ; |pij||2009al ). We have 
computed the uncertainties of the sin i parameters for these 
planetary systems and obtained A(sini) = 0.19, A(sini) = 
3.87 and A (sini) = 6.14 for HD 73526, HD 128311 and 
HD 155358, respectively. The value for HD 73526 well agrees 
with the result of the MCMC simulations, while in the case 
of HD 128311 and HD 155358, these uncertainties are def- 
initely larger than 1, indicating that the amount and/or 
quality of available radial velocity data is not sufficient for 
constraining the orbital inclination. Note that in general, 
Fisher matrix analysis may underestimate the uncertainties 
if the probability distributions of the fitted variables cannot 
be approximated by Gaussian distributions. However, in our 
cases, where the individual measurements are uncorrelated, 
their formal errors are definitely smaller than the amplitude 
of the signal and the sampling of the model function is ade- 
quate (roughly homogeneous for both periods), such a linear 
analysis yields reliable results. 

The results of this analysis have been compared 
with the results provided by Systemic Console package 



|Meschiari et al.ll2009l ) for the planetary system HD 73526. 
By taking into account the mutual perturbations, the two 
applications yielded the same values for the best-fit param- 
eters, however, the residual minimization procedure seemed 
to be more sensitive for the initial parameters in the case of 
the Systemic package. This is mainly due to the inadequate 
choices of the orbital elements 9 The uncertainties derived 
by its built-in bootstrap method were roughly in the same 
magnitude, however, we were unable to derive such a large 
set of points that in the case of If it/nbrv since Systemic 
is slower by roughly two orders of magnitude. 



4.2 Observation scheduling 

Since the program If it is out-of-the-box capable to perform 
analyses on arbitrary user input for which partial derivatives 
are known and have an analytic property (thus, it includes 
the usage of the nbrv_*() functions as well), these features 
can be exploited to optimize observation strategies in or- 
derto derive m ore accu r ate or bital parameters. Recently, 
iFordl (J2008h and iBaluevi (J2008h gives methods with which 
such strategies can be planned efficiently. As it is known 
(see these papers) , the computation of all of the conditional 
probabilities, the expected information content (Ford 2008), 
the D-optimal a nd L-optimal scheduling instances (JBaluevI 
120081 ; IPall l2009al ) requires the evaluation of the covariance 
matrices in arbitrary instances. Since the program If it is 
capable for such an evaluation for arbitrary input functions, 
with the aid of this program, these computations related 
to the optimal strategies can be performed as well without 
any serious difficulties. By default, the program yields both 
the inv e rse of the information matrix Q as defined also by 
IBaluevi (J2008I ) and the goodness sta tistics y 2 , th at and also 
appears in equations (13) - (15) of IFordl (|2008i ). Addition- 
ally, If it is capable to perform these kind of linear analyses 
on arbitrary linear subspace of the parameter domain (i.e. 
parameters in the orthogonal subspace are assumed to be 
fixed or known for independent sources), therefore strategies 
can be built for optimizing various combinations of orbital 
parameters. 



9 The Systemic package employs mean anomaly, eccentricity and 
longitude of pericenter instead of mean longitude and Lagrangian 
elements. A fit performed by If it/nbrv is also more unstable if 
the former set of orbital elements are used. 
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5 SUMMARY 

In this paper we described an algorithm based on the Lie- 
integration method that efficiently computes the paramet- 
ric derivatives of radial velocity model functions for multi- 
ple planetary systems when the planet-planet interactions 
are also taken into account. The analysis of these systems 
yields more accurate constrains for planetary masses since 
the orbital and mutual inclinations can also be derived if 
precise radial velocity data are available. Additionally, the 
presented analytic formulae and integration method aid to 
plan observation schedules in order to optimize the telescope 
time utilization in order to detect planetary perturbations. 
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APPENDIX A: MOTION IN A REFERENCE FRAME FIXED TO ONE OF THE BODIES 

Recalling |Pal fe Siilil (|2007l) . Appendix C, the recurrence relations for the iV-body problem around a fixed center can be 
written as 

L n+1 r irn — L n w im , (Al) 

Li Aijm = Li Tim Li Tjm-, V / 

L n B ijm — L n w im — L n ~w jm , (A3) 

L n K = ^r\L k r lm L n - k w lm , (A4) 

fc=0 V / 
n / \ 
Li Aij = ^/ I . \ -L Aij m L/ ijijm, l-^^J 

fe = V / 

n / \ JV n / \ 

L n+1 «w = -G(X +m,)^ \ \L k <j>iL n ~ k r irn - G £ m^ £ I/ fc <feI/ I_fc Ai.,m + L k 4>jL"~ k r jm , (A6) 

fe=o V / i=i,j^i ft=o V / 

n 

L n+1 <j, t = p- 2 Y, F nkL n - k <t> i L k A i , (A7) 

fc=0 
n 

L" +1 <fe = p^ 2 ^F nk L n -"^ ii L h Aij, (A8) 

where _F n fe = (— 3)(^) + (— 2)( fc ™J. Since the masses (both Al and rrii) are constants, the Lie-derivatives of these are simply 
L n+1 M = L n+1 rrii — for ^ n. Let us denote the linearized of r; m , Wi m , rrii and A4 by £i m , i]im, rrii and 9JT, respectively. 
The vectors H and T> must be extended with the linearized variables m^ and SOT, namely 

S=(«*M*-M«*>,aK) and *=({^}'{4;}'{fl^}'^) (A9) 

and hence 

a '^fe^) + (s^) + (? m 'i) +ffl ^ 

It can easily be shown that the complete set of linearized equations extended with the variables rrii and 97? are 

L n+1 £, lm = L n rn m , (All) 

Li O-ijm = Li C,im — Li qjm, \Al^) 

Li Pijm = Li Tjim Li Tjjm: (AloJ 



E-VL"A Z — y^ I J [L k ^imL n k w im + L k r im L n k rj lm J , (A14) 

E-f>L"Aij = ^2i\lL k aij m L n ~ k Bi jm + L k Aij m L n ~ k l3ij m \, (A15) 

ft=0 V / 

n / \ AT n / \ 

L" +1 7? !m = -G(OT + «i) £ " L k ^L n ~ k r %m -G J2 m ?J2 U [i fe ^ii"" fc ^m + i fe ^L"- fc r jm ] - (A16) 

fc=0 V / j=l,j^ti fc=0 \ / 

-G(M + m.) £ M [(E • VL k fr)L n - k r lm + L*&L n -*fc m ] - 

fc=0 V / 

JV n / \ 

-G 2J nij 2J I j. ] (" " 'DL k (f>ij)L n ~ k Aij m + L k (j} i jL' n ~ k a i j rn + (H ■ T>L k <)>j)L n ~ k rj m + L k (j}jL n ~ k ^ jm , 
S ■ OL" +1 ^ = -2pr 2 ^mr lm L" +1 ^ + ft 2 £ *"»* [(3 ■ VL n ~ k cj>i)L k A t + L n ~ k 'fc(& ■ VL k A t )] , (A17) 

fc=0 

n 

S ■ VL n+1 4> tJ = -2pr 2 Q, jm ^ jm L n+1 ^ + p.^ 2 £ ^nfe [(H ■ ■DL n - k <p ij )L k A ij + L n - k ^j{k ■ t>L k A ij )] . (A18) 

fc=0 

Obviously, the auxiliary variables S]„, S\" m , E^ and S'"^ can be introduced as well fsee lPal fc Sulill2007l . Appendix D), in 
order to optimize the evaluation of equations (|A11|) - ()A18|) . 
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APPENDIX B: PARTIAL DERIVATIVES OF THE COORDINATES AND VELOCITIES 

The computation of equation (|29[) requires the partial derivatives of the initial coordinates and velocities with respect to the 
initial orbital elements. Let us write the initial normalized coordinates and velocities as 



d_(i 



+ 



v 



1 



1 + J 

—s 
+c 



+ 



1 + J 



k 
h 

+h 
-k 



(Bl) 



(B2) 



The normalized coordinates and velocities do not depend on the semimajor axis and the mean motion, therefore these 
quantities are only functions of the mean longitude A and the Lagrangian orbital elements (k,h). In the above equations, 
p = p(A, k,h), q = q(A, k, h), c = cos(p + A), s = sin(p + A) and the quantity J is defined as J = y/1 — e 2 — y/1 — k 2 — h 2 . 

Thus, the partial derivatives of the mass parameter Gm, coordinates (a;, y) and velocities (x, y) with respect to the central 
mass parameter GM and the spectroscopic orbital elements - normalized semi-amplitude K,, mean motion n, mean longitude 
A and the Lagrangian orbital elements (k, h) - are then 



/ 



d(Gm,x,y,x,y) 
d(GM,K,n,X,k,h) 



2 in 



3K. 2 {M+mf 



K. 3 (M+m) 3 



3M + m 

da 



n(3M + m)m 2 n 2 (3M + m)m 2 



8{GM) 

da 
d(GM) 

da 



s 



'/ 



da 

dic^ 

da 



da 
da 



dn 



v 



Hereft'.i/) 

da 
d{GM) 
da 

die 

da 
dn 



dX 



d{GM) 
da 
\d{GM) 

and 



<' 



nij 



da 

die 

da 

die 



,,{' 



I IT) 



da\ . 



da 

a + n— 

dn 




<' 

ail' 

dn' 
m 9A 





l dk 

drj 

l dk 
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The partial derivatives of the normalized coordinates (£,ij) and the normalized velocities (£',»/) with respect to the orbital 
elements (A, k, h) are the following: 
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